co = sns.colors
The study domain and the model domain setup
#
fontsize_tmp = fontsize
fontsize=20
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
topocmap = plt.cm.gist_earth
topocmap.set_bad('w', alpha=0)
#First the greater area map
with nc(os.path.join(geodataf,'topo_2.nc')) as topof:
Gtopo=np.ma.masked_equal(topof.variables['topo'][:],0)
Gtopolon=topof.variables['longitude'][:]
Gtopolat=topof.variables['latitude'][:]
#Now only the tiwi-islands manp
with nc(os.path.join(geodataf,'Topo.nc')) as topof:
topo=np.ma.masked_equal(topof.variables['topo'][:],0)
topolon=topof.variables['longitude'][:]
topolat=topof.variables['latitude'][:]
with nc(CPOLF) as fnc:
lon=fnc.variables['longitude'][:]
lat=fnc.variables['latitude'][:]
fig = plt.figure(figsize=(20,20))
ax = fig.add_subplot(111)
ls = LightSource(azdeg=315, altdeg=45)
M = Basemap(projection='rotpole', llcrnrlat=min(Gtopolat), llcrnrlon=min(Gtopolon), urcrnrlat=max(Gtopolat), urcrnrlon=max(Gtopolon),
resolution='c', area_thresh=1, ax=ax, lon_0=0, o_lon_p=0., o_lat_p=90.,)
borders = (([122.963, 139.055], [-17.091, -7.443]),
([127.385, 134.693], [-14.505, -10.005]),
([129.557, 132.529], [-13.757, -10.757]))
text=(('UM 4km', fontsize+8), ('UM 1.33km', fontsize+8), ('UM 0.44km', fontsize+8))
M.pcolormesh(Gtopolon, Gtopolat, ls.hillshade(Gtopo[:], vert_exag=1), cmap= topocmap, ax=ax)
#M.drawmapscale(M.boundarylons[0]-0.7, M.boundarylats[0]+0.7, max(topolon), min(topolat), 10,
# barstyle='fancy', fontsize=19, labelstyle='simple')
for conf in range(3):
lons,lats = borders[conf]
txt, fts = text[conf]
x, y = M([lons[0], lons[0], lons[1],lons[1]], [lats[0], lats[1], lats[1], lats[0]])
xy = list(zip(x,y))
ttt = plt.text(*M(lons[0]+0.1, lats[0]+0.05), txt, fontsize=fts)
ttt.set_bbox(dict(facecolor='white', alpha=0.7, edgecolor='none', lw=1))
poly = Polygon(xy, alpha=1, edgecolor='k', lw=1, facecolor='none')
plt.gca().add_patch(poly)
#M.drawcoastlines()
M.scatter(*M(131.043101, -12.250323), marker='*',color='k', s=200)
plt.text(*M(borders[0][0][0]+0.15, borders[0][-1][1]-0.6), 'b)', fontsize=fontsize+12)
axins = zoomed_inset_axes(ax, 7, loc=2)
axins.set_xlim(min(topolon), max(topolon))
axins.set_ylim(min(topolat), max(topolat))
tmap = Basemap(projection='cyl',llcrnrlat=min(lat), llcrnrlon=min(lon), urcrnrlat=max(lat), urcrnrlon=max(lon),
resolution='c', area_thresh=1, ax=axins)
tmap.pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap= topocmap)
mark_inset(ax, axins, loc1=3, loc2=4, fc="k", ec='k', lw=0.5)
plt.annotate('a)', xy=(0.01 ,0.92), xycoords='axes fraction', fontsize=fontsize+12)
lons, lats = [min(Gtopolon), max(Gtopolon)], [min(Gtopolat), max(Gtopolat)]
axins2 = inset_axes(ax, width='20%', height='20%', loc=4)
Worldmap= Basemap(projection='ortho', lat_0=Gtopolat[int(len(Gtopolat)/2)], lon_0=Gtopolon[int(len(Gtopolon)/2)],
ax=axins2)
Worldmap.drawcoastlines(color='white', linewidth=0.7)
Worldmap.fillcontinents(color='gray')
bx, by = Worldmap(M.boundarylons, M.boundarylats)
xy = list(zip(bx,by))
poly = Polygon(xy, alpha=1, edgecolor='none', lw=0, facecolor='darkcyan')
_ = Worldmap.ax.add_patch(poly)
_ = fig.subplots_adjust(bottom=0.35)
_ = fig.savefig('Vid/Topo.png', bbox_inches='tight', format='png', dpi=72)
fontsize = fontsize_tmp
Now calculate the map of spectral variance in the diurnal cycle
mpld3.disable_notebook()
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
from matplotlib.ticker import ScalarFormatter
#Plot the output map
fig = plt.figure(figsize=(15,3))
cmap = colmap2.Blues
cmap2 = colmap.GMT_drywet
cmap.set_under('w'),cmap.set_bad('w')
cmap2.set_under('w'),cmap2.set_bad('w')
ax1 = fig.add_subplot(131)
ax1.set_title('Rainfall Climatology', size=fontsize)
im = m1.pcolormesh(lon.filled(-50), lat.filled(-50), map1, vmin=0, cmap=colmap.GMT_drywet, ax=ax1)
m1.drawcoastlines(linewidth=.7)
cbar = m1.colorbar(im , location='bottom', pad="5%")
cbar.set_label('Rain-Rate [mm/d]', size=fontsize)
cbar.ax.tick_params(labelsize=fontsize-4)
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
ax2 = fig.add_subplot(132)
sns.lineplot('x', "y", data=sns_data, ax=ax2, lw=1)
ax2.set(xscale="log")
ax2.set_xlim(0.25,60)
ax2.set_ylim(0,0.6)
ax2.set_xlabel('Periods ($\\tau$) [d]',size=fontsize)
ax2.set_title('Spectral variance [\%]',size=fontsize)
ax2.tick_params(labelsize=fontsize-4)
ax2.set_ylabel('')
ax2.xaxis.set_major_formatter(ScalarFormatter())
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
ax3 = fig.add_subplot(133)
im2 = m1.pcolormesh(lon, lat, map2, vmin=0, vmax=0.8,cmap='Blues', ax=ax3)
m1.drawcoastlines(linewidth=.7)
ax3.set_title('Diurnal Cycle', size=fontsize)
cbar = m1.colorbar(im2,location='bottom',pad="5%")
cbar.set_label('Spectral variance [\%]', size=fontsize)
cbar.ax.tick_params(labelsize=fontsize-4)
# %%
fig.subplots_adjust(right=0.98, bottom=0.03, top=0.99,left=0.01, hspace=0, wspace=0.17)
#'''
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
Burst vs break
# Define extreme Events in the CPOL era
extr1h = pd.read_hdf(extremeTS, 'P10min')
monsoon10m = pd.read_pickle(BurstF)
monsoon10m = monsoon10m.dropna()
monsoon1h = monsoon10m.groupby(pd.TimeGrouper('1h')).mean()
breaks = monsoon10m.loc[monsoon10m == 0].index
bursts = monsoon10m.loc[monsoon10m > 0].index
thresh = perc['10min'][95]
S = pd.Series(np.arange(len(times)),index=times)
print('Bursts: %0.2f %% Breaks: %0.2f %%'
%(len(bursts)/len(monsoon10m.index) * 100, len(breaks)/len(monsoon10m.index) * 100 ))
Bursts: 39.73 % Breaks: 60.27 %
#Plot the diurnal Cycle
#mpld3.enable_notebook()
fig, ax = plt.subplots(1, figsize=(15,10))
fig.subplots_adjust(bottom=0.2, top=0.98, left=0.1, right=0.98)
sns.lineplot(x='Local Time', y='Rain-Rate', data=diurnacl_cycle, hue='Monsoon-Stage', lw=3, ax=ax)
ax.set_ylabel('Rain-Rate mm/d', fontsize=fontsize)
ax.set_xlabel('Local Time [h]', fontsize=fontsize)
#ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%_H'))
ax.legend(loc=0, fontsize=fontsize-4)
ax.tick_params(labelsize=fontsize)
fig.subplots_adjust(bottom=0.145, top=0.99, right=0.99, left=0.1)
The above plots show that diurnal cycle is the most important mode of variability in the considered area. The island thunderstrom hector is the most extreme convective system that occurs at the time-scale at perdiods of 24 h. The Hector 'signal' is most pronounced during the Monsoon break period. Before studying Hector we split the CPOL area into break and burst period and investigate the occurrence of extremes during Monsoon break and burst.
First we define extreme events by applying the 99th percentile threshold of the total 10min dataset. Then we define the break and burst events by the definition of Narsey et al. and calculate the number of burst and the number of break events during the CPOL era
# Now Plot the maps
mpld3.disable_notebook()
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
maxval =[(0,15), (0,15), (-12,12), (0, 150), (0, 150), (-200, 200)]
names=['Breaks', 'Bursts', 'Breaks - Bursts',
'Extremes During Breaks', 'Extremes During Bursts', 'Breaks - Bursts']
cmap_list = [(colmap.GMT_drywet,'k'), (colmap.GMT_drywet,'k'), (colmap.GMT_polar_r,'k'),
(mpl.cm.magma,'w'), (mpl.cm.magma,'w'), (colmap.GMT_polar_r,'k')]
C = [Mean[1]*1.5, Mean[2], Mean[1]*1.5 - Mean[2], D[1]*1.2, D[2], D[1]*1.2 - D[2]]
fig = plt.figure(figsize=(15,10))
mask = ((D[0] * 0 )+1)
im=[]
for i, idx in enumerate(C):
ax = fig.add_subplot(2,3,i+1)
data = np.zeros([len(lat),len(lon)])
colm, coast = cmap_list[i]
colm.set_bad(dict(k='w', w='k')[coast])
colm.set_under(dict(k='w', w='k')[coast])
try:
im.append(m.pcolormesh(lon, lat, mask*C[i].filled(0),vmin=maxval[i][0],vmax=maxval[i][1],cmap=colm))
except AttributeError:
im.append(m.pcolormesh(lon, lat, mask*C[i],vmin=maxval[i][0],vmax=maxval[i][1], cmap=colm))
m.drawcoastlines(color=coast)
ax.set_title('%s'%names[i],size=fontsize)
fig.subplots_adjust(right=0.98, bottom=0.01, top=0.98,left=0.01, hspace=0, wspace=0)
cbar_ax1 = fig.add_axes([0.05, 0.57,0.58, 0.015])
cbar1 = fig.colorbar(im[0], cax=cbar_ax1, orientation='horizontal')
cbar1.ax.tick_params(labelsize=fontsize-4)
cbar1.set_label(u'Rain-Rate [mm/d]',size=fontsize)
cbar_ax2 = fig.add_axes([0.685, 0.57,0.27, 0.015])
cbar2 = fig.colorbar(im[2], cax=cbar_ax2, orientation='horizontal')
cbar2.ax.tick_params(labelsize=fontsize-4)
cbar2.set_label(u'Difference',size=fontsize)
cbar_ax3 = fig.add_axes([0.05, 0.09,0.58, 0.015])
cbar3 = fig.colorbar(im[3], cax=cbar_ax3, orientation='horizontal')
cbar3.ax.tick_params(labelsize=fontsize-4)
cbar3.set_label(u'N. of events above 99th%',size=fontsize)
cbar_ax4 = fig.add_axes([0.685, 0.09,0.27, 0.015])
cbar4 = fig.colorbar(im[-1], cax=cbar_ax4, orientation='horizontal')
cbar4.ax.tick_params(labelsize=fontsize-4)
cbar4.set_ticks([-200, -100, 0, 100, 200])
cbar4.set_label(u'Difference',size=fontsize)
Now we compare the ensemble timeseries of the simulated rainfall for the 1.33km and 0.44km Simulations with the CPOL observations.
time, obs, rain = [], [], []
for e in range(len(ensembles)):
r1 = list(UM133ens.variables['lsrain'][e].mean(axis=(1,2,3)).values)
r2 = list(UM044ens.variables['lsrain'][e].mean(axis=(1,2,3)).values)
obs+=len(r1)*['UM 1.33km'] + len(r2)*['UM 0.44km']
t1 = pd.DatetimeIndex(UM133ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
rain += r1 + r2
t2 = pd.DatetimeIndex(UM044ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
time += list(t1) + list(t2)
r3 = list(OBS.dataset[1].variables['lsrain'].mean(axis=(1,2)).values)
t3 = list(pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime())
o3 = len(r3)*['CPOL']
rain += r3
obs += o3
time += t3
plot_data = pd.DataFrame({'rain-rate':rain, 'time':time, 'Type':obs})
#Plot area avg TS
import matplotlib.dates as mdates
mpld3.disable_notebook()
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
myFmt = mdates.DateFormatter('%d/%m')
#
#Get the minimum time period that is covered by all simulations and also observations
fig = plt.figure(figsize=(15,7))
ax = fig.add_subplot(111)
sns.lineplot(x='time', y='rain-rate', hue='Type', ax=ax, data=plot_data, lw=1.3)
plt.subplots_adjust(left=0.1, bottom=0.1, right=0.98, top=0.9)
ax.legend(loc=0, fontsize=fontsize-4)
ax.tick_params(labelsize=fontsize-4)
ax.set_ylim(-0.001,0.5)
ax.set_ylabel('Rain-Rate [mm/h]', fontsize=fontsize)
_ = ax.xaxis.set_major_formatter(myFmt)
_ = ax.set_xlabel('')
_ = ax.set_title('Comparison of Area-Avg Rainfall', fontsize=fontsize)
fig.subplots_adjust(bottom=0.145, top=0.9, right=0.99, left=0.1)
Now lets look at some animations of UM rainfall occuring over the Tiwi-islands and compare it to the cpol observations
os.chdir('/home/unimelb.edu.au/mbergemann/Work/Programming/Extremes/python/')
#Plot Animatoin
import io
import base64
from IPython.display import HTML
outvid=os.path.join('Vid','WeekOfHector-Ens-3.mp4')
video = io.open(outvid, 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" width="950" height="500" loop="true" controls>
<source src="data:video/mp4;base64,{0}" type="video/mp4" />
</video>'''.format(encoded.decode('ascii')))
#Create the avg. maps of rainfall 0.44km
#m1, m2 = None, None
mpld3.disable_notebook()
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
rom = {1:'i', 2:'ii', 3:'iii', 4:'iv', 5:'v', 6:'vi', 7:'vii', 8:'viii'}
fig = plt.figure(figsize=(15,12))
UM = UM044ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
colm.set_bad('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01) / 6
um_data = np.ma.masked_less(UM.variables['lsrain'][:].values, -0.01) / 6
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM.coords['lat'][:], UM.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
if m1 is None:
m1 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=1)
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m2.pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m2.drawcoastlines(linewidth=0.8)
cbar_ax = fig.add_axes([0.12, 0.12, 0.74, 0.02])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Rain-Rate [mm/d]',size=24)
tit = ['CPOL']
ax[-1].set_title('CPOL', fontsize=fontsize)
for i in range(1,9):
if m2 is None:
m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=1)
ax.append(fig.add_subplot(3,3,i+1))
m2.drawcoastlines(linewidth=0.8)
im.append(m2.pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
tit.append('Ens. %s'%rom[i].upper())
ax[-1].set_title(tit[-1],fontsize=fontsize)
fig.subplots_adjust(right=0.99, bottom=0.15, top=0.95,left=0.01, hspace=0.2, wspace=0.0)
#Create the avg maps of rainfall for 1.33km
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
fig = plt.figure(figsize=(15,12))
UM = UM133ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
colm.set_bad('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01) / 6
um_data = np.ma.masked_less(UM.variables['lsrain'][:].values, -0.01) / 6
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM.coords['lat'][:], UM.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
if m1 is None:
m1 = Basemap(llcrnrlat=min(lat1), llcrnrlon=min(lon1), urcrnrlat=max(lat1), urcrnrlon=max(lon1), resolution='f',
area_thresh=1)
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m2.pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m2.drawcoastlines(linewidth=0.8)
cbar_ax = fig.add_axes([0.12, 0.12, 0.74, 0.02])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Avg. Rain-Rate [mm/d]',size=24)
tit = ['CPOL']
ax[-1].set_title('CPOL', fontsize=fontsize)
for i in range(1,9):
if m2 is None:
m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
area_thresh=1)
ax.append(fig.add_subplot(3,3,i+1))
m2.drawcoastlines(linewidth=0.8)
#m[-1].shadedrelief()
im.append(m2.pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
tit.append('Ens. %s'%rom[i].upper())
ax[-1].set_title(tit[-1],fontsize=fontsize)
fig.subplots_adjust(right=0.99, bottom=0.15, top=0.95,left=0.01, hspace=0.2, wspace=0.0)
#Plet the avg maps for comparison
fig = plt.figure(figsize=(18,12))
#sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
UM_2 = UM133ens.groupby('t.day').sum(dim='t').mean(dim='day')
UM_1 = UM044ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01)
data = [obs_data[1:-1]/6,
np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).mean(axis=0)[0]/6,
np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).mean(axis=0)[0]/6,
None,
np.ma.masked_less(UM_1.variables['lsrain'][:].values/6, -0.01).std(axis=0)[0],
np.ma.masked_less(UM_2.variables['lsrain'][:].values/6, -0.01).std(axis=0)[0]
]
data.append(None)
data.append(data[0] - data[1])
data.append(data[0] - data[2])
UM_1min = np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).min(axis=0)[0]
UM_2min = np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).min(axis=0)[0]
UM_1max = np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).max(axis=0)[0]
UM_2max = np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).max(axis=0)[0]
#data.append(None)
#data.append(UM_1max - UM_1min)
#data.append(UM_1max - UM_1min)
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM_1.coords['lat'][:], UM_1.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m2.pcolormesh(lon1, lat1, obs_data/6,vmin=0.0,vmax=25/6,cmap=colm)]
m2.drawcoastlines()
#cbar_ax = fig.add_axes([0.14, 0.33, 0.74, 0.01])
#cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
#cbar.ax.tick_params(labelsize=24)
#cbar.set_label('Avg. Rain-rate [mm/d]',size=24)
tit = ['CPOL avg.', 'UM 0.44km ens avg.', 'UM 1.33km ens avg.',
None, 'UM 0.44 ens std.', 'UM 1.33km ens std.',
None, 'CPOL - UM 0.44 ens avg.' ,'CPOL - UM 1.33km ens avg.',
None, 'UM 0.44 ens range', 'UM 1.33km ens range']
colors = {1:colm, 2:colm, 4:colm, 5:colm,
7:colmap.GMT_polar_r, 8:colmap.GMT_polar_r,
10:colm, 11:colm}
Range = {1:(0.0, 25/6), 2:(0.0, 25/6), 4:(0,8/6), 5:(0,8/6), 7:(-20/6,20/6), 8:(-20/6,20/6), 10:(0,30/6), 11:(0,30/6)}
height = {2:0.682, 5:0.37, 8:0.044, 11:0.05}
ax[-1].set_title('CPOL avg.', fontsize=fontsize-2)
cbar_ax = [fig.add_axes([0.91, height[2], 0.01, 0.23])]
cbar = [fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical')]
cbar[-1].ax.tick_params(labelsize=fontsize-4)
cbar[-1].set_label('Rain-Rate [mm/d]',size=fontsize-4)
for i in range(1,len(data)):
if data[i] is not None:
ax.append(fig.add_subplot(3,3,i+1))
#m[-1].shadedrelief()
im.append(m2.pcolormesh(lon2, lat2, data[i],vmin=Range[i][0],vmax=Range[i][1],cmap=colors[i]))
m2.drawcoastlines()
ax[-1].set_title(tit[i],fontsize=fontsize-2)
if i in (5, 8, 11):
cbar_ax.append(fig.add_axes([0.91, height[i], 0.01, 0.23]))
cbar.append(fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical'))
cbar[-1].ax.tick_params(labelsize=fontsize-4)
cbar[-1].set_label('Rain-rate [mm/d]',size=fontsize-2)
fig.subplots_adjust(right=0.9, bottom=0.01, top=0.95,left=0.01, hspace=0.0, wspace=0)
data, time, model = [], [], []
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
for i in range(len(ensembles)):
d1 = UM133.dataset[i].groupby('t.hour').mean(dim='t').mean(dim=('surface','lat','lon'))
d2 = UM044.dataset[i].groupby('t.hour').mean(dim='t').mean(dim=('surface','lat','lon'))
t1 = pd.Series(d1.variables['lsrain'].values, index=(d1.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
t2 = pd.Series(d2.variables['lsrain'].values, index=(d2.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
time += list(t1.index) + list(t2.index)
data += list(t1.values) + list(t2.values)
model += len(t1)*['UM 1.33km'] + len(t2)*['UM 0.44km']
cpol_1 = OBS.dataset[1].groupby('t.hour').mean(dim='t').mean(dim=('y','x'))
obs_S = pd.Series(cpol_1.variables['lsrain'].values, index=(cpol_1.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
data += list(obs_S.values)
time += list(obs_S.index)
model += len(obs_S)*['Cpol']
plot_data = pd.DataFrame(dict(data=data, time=time, Type=model))
#mpld3.enable_notebook()
# Plot avg diurnal cycle
#Get the minimum time period that is covered by all simulations and also observations
fig = plt.figure(figsize=(15,8), dpi=72)
ax = fig.add_subplot(111)
####################
sns.lineplot(x=time, y=data, hue=model, ax=ax, lw=3, data=plot_data)
#ax.set_ylim(0.0,0.42)
ax.set_xlim(0,23)
ax.set_ylabel('Rain-Rate [mm/h]', fontsize=fontsize)
ax.set_xlabel('Local Time [h]', fontsize=fontsize)
ax.legend(loc=0, fontsize=fontsize-4)
ax.tick_params(labelsize=fontsize-4)
fig.subplots_adjust(bottom=0.145, top=0.99, right=0.99, left=0.1)
#ax.grid(color='r', linestyle='-', linewidth=0)
mpld3.disable_notebook()
#Plot the maps
plt.close()
mpld3.disable_notebook()
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':False})
um_1 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('surface'))
um_2 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('surface'))
obs = OBS.dataset[1].groupby('t.hour').mean(dim='t')
time = (um_1.hour + 9) % 24
data = [obs, um_2, um_1]
fig = plt.figure(figsize=(15,7))
lat1, lon1 = OBS.dataset[1].variables['latitude'].values[:,0], OBS.dataset[1].variables['longitude'].values[0,:]
lat2, lon2 = um_1.coords['lat'][:], um_1.coords['lon'][:]
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.98,left=0.01, hspace=0, wspace=0)
first = True
val_max = 2
colm = colmap2.Blues
colm.set_under('w', alpha=0)
colm.set_bad('w', alpha=0)
#ax = [fig.add_subplot(2,3,1)]
mm, im = [], []
col = 2
for i in range(24):
#break
fname_1 = os.path.join(outdir, 'DiurnalCycle_%02i.png'%i)
if first:
for col in range(1, 7):
if col != 4:
ax1 = fig.add_subplot(2,3,col)
m2.drawcoastlines(linewidth=0.8)
if col == 1:
ax1.set_title('CPOL',fontsize=fontsize)
_ = im.append(m2.pcolormesh(lon1, lat1, obs['lsrain'].values[i], shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
elif col == 2:
ax1.set_title('UM 1.33km mean', fontsize=fontsize)
_ = im.append(m2.pcolormesh(lon2, lat2, um_2['lsrain'].mean(dim='ens').values[i], shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
elif col == 3:
ax1.set_title('UM 0.44km mean', fontsize=fontsize)
_ = im.append(m2.pcolormesh(lon2, lat2, um_1['lsrain'].mean(dim='ens').values[i], shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
elif col == 5:
ax1.set_title('UM 1.33km std', fontsize=fontsize)
_ = im.append(m2.pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i], shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
elif col == 6:
ax1.set_title('UM 0.44km std', fontsize=fontsize)
_ = im.append(m2.pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i], shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
cbar_ax = fig.add_axes([0.14, 0.0, 0.74, 0.02])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
_ = cbar.ax.tick_params(labelsize=fontsize-4)
first = False
else:
###MEAN
im[0].set_array(obs['lsrain'].values[i].ravel())
im[1].set_array(um_2['lsrain'].mean(dim='ens').values[i].ravel())
im[2].set_array(um_1['lsrain'].mean(dim='ens').values[i].ravel())
###STD
im[-2].set_array(um_2['lsrain'].std(dim='ens').values[i].ravel())
im[-1].set_array(um_1['lsrain'].std(dim='ens').values[i].ravel())
cbar.set_label('Rain-Rate at %2i LT [mm/h]'%((i+9.5)%24),size=fontsize)
fig.savefig(fname_1, bbox_inches='tight', format='png', dpi=72)
#break
dest_dir = os.path.abspath('Vid')
make_mp4_from_frames(outdir, dest_dir, 'WeekOfHector-Diurnal-2.mp4', 4, glob='DiurnalCycle_??')
fig.clf()
plt.close()
day = (time[1:5] - 9) % 24
night = (time[5:5+4] - 9) % 24
obs_diff = obs['lsrain'][day].sum(dim='hour') - obs['lsrain'][night].sum('hour')
um_1_diff = um_1['lsrain'][day].sum(dim='hour') - um_1['lsrain'][night].sum('hour')
um_2_diff = um_2['lsrain'][day].sum(dim='hour') - um_2['lsrain'][night].sum('hour')
mask = np.ma.masked_invalid(obs['lsrain'].values) * 0 + 1
obs_max = (((np.argmax(obs['lsrain'].values * mask, axis=0)*mask[0] + 9.5 ) % 24)/2).round(0) * 2
um_1_max = ((um_1['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9.5) % 24 / 2).round(0) * 2
um_2_max = ((um_2['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9.5) % 24 / 2).round(0) * 2
colm = CubeYF_20.get_mpl_colormap(N=24, gamma=2.0)
#colm = qualitative.Paired[12].get_mpl_colormap(N=24, gamma=2.0)
colm.set_bad('w')
#m, ax = [], []
ax = []
data = [obs_max, um_1_max, um_2_max]
names = ['CPOL', 'UM 1.33km', 'UM 0.44km']
fig = plt.figure(figsize=(15,5.4))
fig.subplots_adjust(right=0.98, bottom=0.03, top=0.98,left=0.01, hspace=0, wspace=0)
cbar_ax = fig.add_axes([0.125, 0.15, 0.75, 0.05])
for i in range(3):
ax.append(fig.add_subplot(1,3,i+1))
if m2 is None:
m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
area_thresh=1)
m2.drawcoastlines()
try:
im = m2.pcolormesh(lon1, lat1, data[i], vmin=0, vmax=24, cmap=colm)
except TypeError:
im = m2.pcolormesh(lon2, lat2, data[i], vmin=0, vmax=24, cmap=colm)
ax[-1].set_title(names[i], fontsize=fontsize)
cbar=fig.colorbar(im, cax=cbar_ax, orientation='horizontal',ticks=list(range(1,24,2)))
cbar.ax.tick_params(labelsize=fontsize-4)
cbar.set_label('Local Time [h]', size=fontsize)
mpld3.disable_notebook()
#Plot the pdfs
sns.set_style('darkgrid')
from mpl_toolkits.axes_grid.inset_locator import inset_axes
dist=lognorm([LogNorm.CPOL['shape']*3.4], loc = np.log(LogNorm.CPOL['scale']*1.9))#, shape # mu, sigma
med = np.median(um044vec.loc[um044vec>0].values)
X = np.linspace(0,65,200)
#b = 2.5
nbins = 100
def fitpdf(x, b):
a=med*100
return ((b/a) * (x/a)**(b-1)) / (1+(x/a)**b)**2
#popt, pcov = curve_fit(powerlaw, um044cnt['Rain-rate'].values[1:],
# um044cnt['count'].values[1:])
popt, pcov = PowerFit.CPOL['popt']/4, PowerFit.CPOL['pcov']
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111)
ax2 = ax.twinx()
#Define histogram data to plot
histdata = (um133vec.loc[um133vec>0], um044vec.loc[um044vec>0], obs_vec.loc[obs_vec>0])
#facecolors = ['fuchsia', 'firebrick', 'lightseagreen']
#n1_1, bins1_1, patches1_1 = ax2.hist(, bins=nbins, normed=1, facecolor='firebrick', alpha=0.25)
#n2, bins2, patches2 = ax2.hist(, bins=nbins, normed=1, facecolor='lightseagreen', alpha=0.25)
ax2.set_ylim(0,0.9)
ax2.set_yticks([])
#ax.plot(np.linspace(0.01, 100,200), fitpdf(np.linspace(0.01,100,200),popt[0]), lw=4, color='navy',
# label='Log-logistic $\\beta = %02.2e$'%popt[0])
#ax.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
#ax.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'--', label='UM 1.33 km ens',lw=4, color='fuchsia')
#ax.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=4, color='lightseagreen')
#ax.scatter(um044cnt['Rain-rate'].values, um044cnt['count'].values, c='red', s=125, alpha=0.24)
#ax.scatter(obscnt['Rain-rate'].values, obscnt['count'].values, c='purple', s=125, alpha=0.24)
#ax.plot(np.linspace(0, 50,200), um044pdf, lw=4, label='UM 0.44km ens', color='lightseagreen')
#ax.plot(np.linspace(0, 50,200), obspdf, lw=4, label='CPOL', color='fuchsia')
#ax.plot(np.linspace(0.0,50,200), dist.pdf(np.linspace(0.0, 50,200)), color='firebrick', alpha=0.8,
# lw=4, label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
n1, bins1, patches1 = ax.hist(histdata, bins=nbins, normed=1, alpha=0.55)
#fitted_pl.plot_pdf(ax=ax,color='k',lw=4)
#ax.plot(np.linspace(0.1, 100), um044dens_aavg(np.linspace(0.1,100)), lw=4, label='UM 0.44km area-avg')
#ax.plot(np.linspace(0.1, 100), um133dens_aavg(np.linspace(0.1,100)), lw=4, label='UM 1.33km area-avg')
ax.tick_params(labelsize=fontsize-4)
#ax.set_ylim(0.,.25)
ax.set_xlim(0.,20)
ax.set_ylabel('Density', fontsize=fontsize)
ax.set_xlabel('Rain-Rate [mm/h]', fontsize=fontsize)
axin = inset_axes(ax, width = "80%", height = "80%", loc=1)
axin.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'-', label='UM 1.33 km ens',lw=2)
axin.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=2)
axin.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=2)
#axin.plot(np.linspace(0.1,80,200), dist.pdf(np.linspace(0.1, 80,200)), lw=1, alpha=0.8,
# label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
#axin.set_xlabel('$\\log$(Rain-rate)', fontsize=fontsize-2)
#axin.plot(X, um044pdf,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
#axin.plot(X, obspdf,'-', label='CPOL',lw=4, color='lightseagreen')
axin.set_yscale('log')
axin.set_xscale('log')
#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=fontsize-4)
axin.set_xlim(0.1, 100)
axin.set_ylim(0.000023, 1)
ann = axin.annotate("(b)", (0.95,0.93), xycoords="axes fraction", fontsize=fontsize)
axin.legend(loc=3, fontsize=fontsize-4)
#n.set_ylabel('$\\log$(Density)', fontsize=fontsize-2)
fig.subplots_adjust(bottom=0.13, right=0.98, top=0.97)
#mpld3.enable_notebook()
Now get the percentiles for the ensemble simulation and compare it against the observations
#Plot the percentages
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111)
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
ax.plot(PERC.index, PERC['UM 1.33km'].values, label='UM 1.33km',lw=3)
ax.plot(PERC.index, PERC['UM 0.44km'].values, label='UM 0.44km', lw=3)
ax.plot(PERC.index, PERC['Obs'].values, label='CPOL',lw=2)
ax.set_xlim(1,100)
#ax.set_ylim(10,100)
#ax.set_yscale('log')
#ax.set_xscale('log')
ax.set_xlabel('Percentile [\%]', fontsize=fontsize)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=fontsize)
ax.tick_params(labelsize=fontsize-4)
axin = inset_axes(ax, width = "80%", height = "60%", loc=9)
axin.plot(PERC.index[1:], PERC['UM 1.33km'].values[1:], label='UM 1.33km',lw=3)
axin.plot(PERC.index[1:], PERC['UM 0.44km'].values[1:], label='UM 0.44km', lw=3)
axin.plot(PERC.index[1:], PERC['Obs'].values[1:], label='CPOL',lw=2)
#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=fontsize-4)
axin.set_xlim(50, 100)
axin.set_ylim(.25,100)
#axin.set_xticks([20, 200, 500])
axin.legend(loc=0, fontsize=fontsize-4)
#axin.get_xaxis().set_major_formatter(mpl.ticker.ScalarFormatter(useMathText=False))
#axin.get_xaxis().get_major_formatter().labelOnlyBase = True
#axin.set_xticklabels([], rotation=0)
ann = axin.annotate("(b)", (0.95,0.93), xycoords="axes fraction", fontsize=fontsize)
#axin.set_xticks(range(50,100,10), minor=False)
axin.set_yscale('log')
axin.set_xscale('log')
fig.subplots_adjust(bottom=0.13, right=0.98, top=0.97)
mpl.ticker.ScalarFormatter?